# Working Directory
setwd("~/Desktop/PolicingTStops/BadApples/Replication/")

# Libraries
library(readr)
library(dplyr)
library(ggplot2)

# Loading in the data.
load("Data/CT_WideOfficerData.RData")
load("Data/CT_FullFile.RData")

# Identify high disparity officers. 
ct.officer.wide$ratio = abs(ct.officer.wide$fasttt.mean_White/ct.officer.wide$fasttt.mean_Black)
summary(ct.officer.wide$ratio)
ct.officer.extreme.bw = boxplot(ct.officer.wide$ratio)
ct.officer.extreme.id = ct.officer.wide$officer[which(ct.officer.wide$ratio%in%
                                                        unique(ct.officer.extreme.bw$out))]
ct$exclude.officer = ifelse(ct$unique.officer %in% ct.officer.extreme.id, 1, 0)
table(ct$exclude.officer)

ct.agencies = unique(apply(as.matrix(as.character(ct.officer.wide$officer)),
                           1,
                           function(x){strsplit(x,"::",fixed=T)[[1]][1]}))
ct.agencies.extreme = unique(apply(as.matrix(as.character(ct.officer.wide$officer[which(ct.officer.wide$ratio%in%unique(ct.officer.extreme.bw$out))])),
                           1,
                           function(x){strsplit(x,"::",fixed=T)[[1]][1]}))

# Collapsing by Agency.
ct.dept.temp = aggregate(ct[,c("stop.white","search.white","contra.white",
                               "stop.black","search.black","contra.black")],
                         by=list(ct$`Department Name`),sum,na.rm=T)
ct.dept = aggregate(ct[,c("stopoccur","searchoccur","contraband")],
                    by=list(ct$`Department Name`,ct$race.ethnicity),
                    sum,na.rm=T)
ct.dept.ids = ct.dept.temp$Group.1[ct.dept.temp$stop.white>49&
                                     ct.dept.temp$stop.black>49&
                                     ct.dept.temp$search.white>0&
                                     ct.dept.temp$search.black>0]
ct.dept$include = ifelse(ct.dept$Group.1%in%ct.dept.ids,1,0)
ct.dept.ids2 = ct.dept.temp$Group.1[ct.dept.temp$stop.white>49&
                                      ct.dept.temp$stop.black>49]
ct.dept$include2 = ifelse(ct.dept$Group.1%in%ct.dept.ids2,1,0)

save(ct.dept,file="Data/CT_Depts_50.RData")

# Collapsing by Agency w/out "extremes."
ct.dept.temp = aggregate(ct[ct$exclude.officer==0,
                            c("stop.white","search.white","contra.white",
                               "stop.black","search.black","contra.black")],
                         by=list(ct$`Department Name`[ct$exclude.officer==0]),
                         sum,na.rm=T)
ct.dept = aggregate(ct[ct$exclude.officer==0,
                       c("stopoccur","searchoccur","contraband")],
                    by=list(ct$`Department Name`[ct$exclude.officer==0],
                            ct$race.ethnicity[ct$exclude.officer==0]),
                    sum,na.rm=T)
ct.dept.ids = ct.dept.temp$Group.1[ct.dept.temp$stop.white>49&
                                     ct.dept.temp$stop.black>49&
                                     ct.dept.temp$search.white>0&
                                     ct.dept.temp$search.black>0]
ct.dept$include = ifelse(ct.dept$Group.1%in%ct.dept.ids,1,0)
ct.dept.ids2 = ct.dept.temp$Group.1[ct.dept.temp$stop.white>49&
                                      ct.dept.temp$stop.black>49]
ct.dept$include2 = ifelse(ct.dept$Group.1%in%ct.dept.ids2,1,0)

save(ct.dept,file="Data/CT_Depts_50_NoOutliers.RData")

##
## Agency Models -- Full
##

rm(list=ls())

# Load in the collapsed data sets. 
load("Data/CT_Depts_50.RData")

# Setting up the officer stan list and cropping the data set further
ct.dept.sm = subset(ct.dept, ct.dept$include==1)
table(ct.dept.sm$include,ct.dept.sm$include2)

ct.dept.sm$state.agency = ifelse(grepl("CSP ",ct.dept.sm$Group.1),1,0)
ct.dept.sm$state.agency = ifelse(grepl("DMV ",ct.dept.sm$Group.1),1,
                                    ct.dept.sm$state.agency)
ct.dept.sm$state.agency = ifelse(grepl("DMV",ct.dept.sm$Group.1),1,
                                    ct.dept.sm$state.agency)
ct.dept.sm$state.agency = ifelse(grepl("MTA",ct.dept.sm$Group.1),1,
                                    ct.dept.sm$state.agency)
ct.dept.sm$state.agency = ifelse(grepl("Yale",ct.dept.sm$Group.1),1,
                                    ct.dept.sm$state.agency)
ct.dept.sm$state.agency = ifelse(grepl("State Police",ct.dept.sm$Group.1),1,
                                    ct.dept.sm$state.agency)
ct.dept.sm$state.agency = ifelse(grepl("SU",ct.dept.sm$Group.1),1,
                                    ct.dept.sm$state.agency)
ct.dept.sm$state.agency = ifelse(grepl("UCONN",ct.dept.sm$Group.1),1,
                                    ct.dept.sm$state.agency)
ct.dept.sm$state.agency = ifelse(grepl("CAPITOL POLICE",ct.dept.sm$Group.1),1,
                                    ct.dept.sm$state.agency)
ct.dept.sm$state.agency = ifelse(grepl("Tribal",ct.dept.sm$Group.1),1,
                                    ct.dept.sm$state.agency)

table(ct.dept.sm$state.agency)
ct.dept.sm = subset(ct.dept.sm,ct.dept.sm$state.agency == 0)
table(ct.dept.sm$state.agency)

ct.dept.sm = subset(ct.dept.sm, ct.dept.sm$Group.2=="White"|
                         ct.dept.sm$Group.2=="Black")

table(ct.dept.sm$include,ct.dept.sm$include2)

colnames(ct.dept.sm) = c("agency","driver_rg","stops","searches","contraband",
                            "include","include2","state.agency")

save(ct.dept.sm, file="Data/CT_AgencySm_Data_50Thresh.RData")

# Agency set up for running Stan model. 
ct.agency.stan = list(N = nrow(ct.dept.sm),
                       R = length(unique(ct.dept.sm$driver_rg)),
                       D = length(unique(ct.dept.sm$agency)),
                       r = as.integer(as.factor(ct.dept.sm$driver_rg)),
                       d = as.integer(as.factor(ct.dept.sm$agency)),
                       n = ct.dept.sm$stops,
                       s = ct.dept.sm$searches,
                       h = ct.dept.sm$contraband)

# The original specification is 5 chains, a 4000 warmup iteaations, 
# 12,000 iterations that are kept, and a delta of 0.95.
library(rstan)
rstan_options(auto_write = TRUE)

options(mc.cores = 5)
ct.model.multi = stan("Scripts/fasttt-master/stan_models/model_mixture.stan",
                      data = ct.agency.stan,
                      chains = 5,
                      cores = 5,
                      warmup = 5000,#Orig: 4000
                      iter = 25000,#Oirg: 12500
                      thin = 5,
                      seed = 98713,
                      control = list(adapt_delta = 0.9,#Orig: 0.95
                                     max_treedepth = 20))

save(ct.model.multi,file="Models/CT_FastTT_Agency50Thresh.RData")


##
## Agency Models -- No Outliers
##

rm(list=ls())

# Load in the collapsed data sets. 
load("Data/CT_Depts_50_NoOutliers.RData")

# Setting up the officer stan list and cropping the data set further
ct.dept.sm = subset(ct.dept, ct.dept$include==1)
table(ct.dept.sm$include,ct.dept.sm$include2)

ct.dept.sm$state.agency = ifelse(grepl("CSP ",ct.dept.sm$Group.1),1,0)
ct.dept.sm$state.agency = ifelse(grepl("DMV ",ct.dept.sm$Group.1),1,
                                 ct.dept.sm$state.agency)
ct.dept.sm$state.agency = ifelse(grepl("DMV",ct.dept.sm$Group.1),1,
                                 ct.dept.sm$state.agency)
ct.dept.sm$state.agency = ifelse(grepl("MTA",ct.dept.sm$Group.1),1,
                                 ct.dept.sm$state.agency)
ct.dept.sm$state.agency = ifelse(grepl("Yale",ct.dept.sm$Group.1),1,
                                 ct.dept.sm$state.agency)
ct.dept.sm$state.agency = ifelse(grepl("State Police",ct.dept.sm$Group.1),1,
                                 ct.dept.sm$state.agency)
ct.dept.sm$state.agency = ifelse(grepl("SU",ct.dept.sm$Group.1),1,
                                 ct.dept.sm$state.agency)
ct.dept.sm$state.agency = ifelse(grepl("UCONN",ct.dept.sm$Group.1),1,
                                 ct.dept.sm$state.agency)
ct.dept.sm$state.agency = ifelse(grepl("CAPITOL POLICE",ct.dept.sm$Group.1),1,
                                 ct.dept.sm$state.agency)
ct.dept.sm$state.agency = ifelse(grepl("Tribal",ct.dept.sm$Group.1),1,
                                 ct.dept.sm$state.agency)

table(ct.dept.sm$state.agency)
ct.dept.sm = subset(ct.dept.sm,ct.dept.sm$state.agency == 0)
table(ct.dept.sm$state.agency)

ct.dept.sm = subset(ct.dept.sm, ct.dept.sm$Group.2=="White"|
                      ct.dept.sm$Group.2=="Black")

table(ct.dept.sm$include,ct.dept.sm$include2)

colnames(ct.dept.sm) = c("agency","driver_rg","stops","searches","contraband",
                         "include","include2","state.agency")

save(ct.dept.sm, file="Data/CT_AgencySm_Data_50Thresh_NoOut.RData")

# Agency set up for running Stan model. 
ct.agency.stan = list(N = nrow(ct.dept.sm),
                      R = length(unique(ct.dept.sm$driver_rg)),
                      D = length(unique(ct.dept.sm$agency)),
                      r = as.integer(as.factor(ct.dept.sm$driver_rg)),
                      d = as.integer(as.factor(ct.dept.sm$agency)),
                      n = ct.dept.sm$stops,
                      s = ct.dept.sm$searches,
                      h = ct.dept.sm$contraband)

# The original specification is 5 chains, a 4000 warmup iteaations, 
# 12,000 iterations that are kept, and a delta of 0.95.
library(rstan)
rstan_options(auto_write = TRUE)

options(mc.cores = 5)
ct.model.multi = stan("Scripts/fasttt-master/stan_models/model_mixture.stan",
                      data = ct.agency.stan,
                      chains = 5,
                      cores = 5,
                      warmup = 5000,#Orig: 4000
                      iter = 25000,#Oirg: 12500
                      thin = 5,
                      seed = 98713,
                      control = list(adapt_delta = 0.9,#Orig: 0.95
                                     max_treedepth = 20))

save(ct.model.multi,file="Models/CT_FastTT_Agency50Thresh_NoOut.RData")

###
### Analysis
###

rm(list=ls())

#
# Building the baseline data set
#

# Load in the collapsed data sets. 
load("Data/CT_AgencySm_Data_50Thresh.RData")
load("Models/CT_FastTT_Agency50Thresh.RData")

# Traceplot for a random sample of thresholds. 
set.seed(987)
t_for_trace = sample(1:176,8)
t_names_for_trace = paste("t_i[",t_for_trace,"]",sep="")

png("Figures/CT_Traceplot_AgencyFull50.png",1254,769)
traceplot(ct.model.multi,
          pars=t_names_for_trace,nrow=2)
dev.off()

png("Figures/CT_SampOffDensity_AgencyFull50.png",629,592)
plot(ct.model.multi, show_density=TRUE, 
     pars=t_names_for_trace, 
     ci_level=0.95, outer_level=0.999)
dev.off()

# Extracting the information + final check.
ct.summary = summary(ct.model.multi)$summary

ct.summary[which(is.na(ct.summary[,10])),]
table(ct.summary[,10]<1)
table(ct.summary[,10]<1.01)

# Cropping the model output to what is needed.
ct.summary.ti = ct.summary[2:177,]
ct.output = data.frame("fasttt.mean" = ct.summary.ti[,1],
                       "fasttt.lower" = ct.summary.ti[,1]-
                         1.96*ct.summary.ti[,2],
                       "fasttt.upper" = ct.summary.ti[,1]-
                         1.96*ct.summary.ti[,2])

# Merging in that information. 
ct.agency.analysis = cbind(ct.dept.sm,ct.output)

# Transforming the data set from long to wide. 
ct.agency.sm.temp = ct.agency.analysis[,c(1:5,9:11)]
ct.agency.wide = tidyr::pivot_wider(data=ct.agency.sm.temp,
                                     id_cols="agency",
                                     names_from="driver_rg",
                                     values_from = c("stops","searches","contraband",
                                                     "fasttt.mean","fasttt.lower",
                                                     "fasttt.upper"))

# Calculating additional values of interest.
ct.agency.wide$search_rate_white = ct.agency.wide$searches_White/
  ct.agency.wide$stops_White
ct.agency.wide$search_rate_black = ct.agency.wide$searches_Black/
  ct.agency.wide$stops_Black
ct.agency.wide$hit_rate_white = ifelse(ct.agency.wide$searches_White>0&
                                          ct.agency.wide$searches_Black>0,
                                        ct.agency.wide$contraband_White/
                                          ct.agency.wide$searches_White,NA)
ct.agency.wide$hit_rate_black = ifelse(ct.agency.wide$searches_Black>0&
                                          ct.agency.wide$searches_White>0,
                                        ct.agency.wide$contraband_Black/
                                          ct.agency.wide$searches_Black,NA)
ct.agency.wide$dif.tt = ct.agency.wide$fasttt.mean_White-
  ct.agency.wide$fasttt.mean_Black
ct.agency.wide$dif.tt.cat = ifelse(ct.agency.wide$dif.tt>0,
                                    ifelse(ct.agency.wide$fasttt.lower_White>
                                             ct.agency.wide$fasttt.upper_Black,
                                           "3 Black Drivers Disp. Neg.",
                                           "2 Approx. Equal"),
                                    ifelse(ct.agency.wide$fasttt.lower_Black>
                                             ct.agency.wide$fasttt.upper_White,
                                           "1 White Drivers Disp. Neg.",
                                           "2 Approx. Equal"))
ct.agency.wide$dif.sr = ct.agency.wide$search_rate_white-
  ct.agency.wide$search_rate_black
ct.agency.wide$dif.sr.cat = ifelse(ct.agency.wide$dif.sr==0,
                                    "2 Apprrox. Equal",
                                    ifelse(ct.agency.wide$dif.sr>0,
                                           "1 White Drivers Disp. Neg.",
                                           "3 Black Drivers Disp. Neg."))
ct.agency.wide$dif.hr = ct.agency.wide$hit_rate_white-
  ct.agency.wide$hit_rate_black
ct.agency.wide$dif.hr.cat = ifelse(ct.agency.wide$dif.hr==0,
                                    "2 Approx. Equal",
                                    ifelse(ct.agency.wide$dif.hr<0,
                                           "1 White Drivers Disp. Neg.",
                                           "3 Black Drivers Disp. Neg."))
colnames(ct.agency.wide)[2:23] = paste(colnames(ct.agency.wide)[2:23],"full",sep="_")

save(ct.agency.wide,file="Data/CT_WideAgencyData_50Full.RData")

ct.threshold.plot = ggplot(ct.agency.wide,aes(fasttt.mean_White_full,fasttt.mean_Black_full)) + 
  geom_point(size=1.5,shape=1) + theme_bw(base_size=12) +
  xlab(NULL) + ylab(NULL) +
  geom_abline(intercept = 0, slope = 1,color = "black",size = 0.5) +
  geom_abline(intercept = 0, slope = 2,color = "darkgray",linetype="dashed",size = 0.5) +
  geom_abline(intercept = 0, slope = 0.5,color = "darkgray",linetype="dashed",size = 0.5) +
  xlim(0,5)+ylim(0,5)+
  theme(panel.grid.major = element_line(colour = NA),
        panel.grid.minor = element_line(colour=NA),
        legend.position = "none") +
  ggtitle("Inferred Thresholds")
ct.threshold.plot

#
# Building the baseline data set
#

# Load in the collapsed data sets. 
load("Data/CT_AgencySm_Data_50Thresh_NoOut.RData")
load("Models/CT_FastTT_Agency50Thresh_NoOut.RData")

# Traceplot for a random sample of thresholds. 
set.seed(987)
t_for_trace = sample(1:176,8)
t_names_for_trace = paste("t_i[",t_for_trace,"]",sep="")

png("Figures/CT_Traceplot_AgencyNoOut50.png",1254,769)
traceplot(ct.model.multi,
          pars=t_names_for_trace,nrow=2)
dev.off()

png("Figures/CT_SampOffDensity_AgencyNoOut50.png",629,592)
plot(ct.model.multi, show_density=TRUE, 
     pars=t_names_for_trace, 
     ci_level=0.95, outer_level=0.999)
dev.off()

# Extracting the information + final check.
ct.summary = summary(ct.model.multi)$summary

ct.summary[which(is.na(ct.summary[,10])),]
table(ct.summary[,10]<1)
table(ct.summary[,10]<1.01)

# Cropping the model output to what is needed.
ct.summary.ti = ct.summary[2:177,]
ct.output = data.frame("fasttt.mean" = ct.summary.ti[,1],
                       "fasttt.lower" = ct.summary.ti[,1]-
                         1.96*ct.summary.ti[,2],
                       "fasttt.upper" = ct.summary.ti[,1]-
                         1.96*ct.summary.ti[,2])

# Merging in that information. 
ct.agency.analysis = cbind(ct.dept.sm,ct.output)

# Transforming the data set from long to wide. 
ct.agency.sm.temp = ct.agency.analysis[,c(1:5,9:11)]
ct.agency.wide = tidyr::pivot_wider(data=ct.agency.sm.temp,
                                    id_cols="agency",
                                    names_from="driver_rg",
                                    values_from = c("stops","searches","contraband",
                                                    "fasttt.mean","fasttt.lower",
                                                    "fasttt.upper"))

# Calculating additional values of interest.
ct.agency.wide$search_rate_white = ct.agency.wide$searches_White/
  ct.agency.wide$stops_White
ct.agency.wide$search_rate_black = ct.agency.wide$searches_Black/
  ct.agency.wide$stops_Black
ct.agency.wide$hit_rate_white = ifelse(ct.agency.wide$searches_White>0&
                                         ct.agency.wide$searches_Black>0,
                                       ct.agency.wide$contraband_White/
                                         ct.agency.wide$searches_White,NA)
ct.agency.wide$hit_rate_black = ifelse(ct.agency.wide$searches_Black>0&
                                         ct.agency.wide$searches_White>0,
                                       ct.agency.wide$contraband_Black/
                                         ct.agency.wide$searches_Black,NA)
ct.agency.wide$dif.tt = ct.agency.wide$fasttt.mean_White-
  ct.agency.wide$fasttt.mean_Black
ct.agency.wide$dif.tt.cat = ifelse(ct.agency.wide$dif.tt>0,
                                   ifelse(ct.agency.wide$fasttt.lower_White>
                                            ct.agency.wide$fasttt.upper_Black,
                                          "3 Black Drivers Disp. Neg.",
                                          "2 Approx. Equal"),
                                   ifelse(ct.agency.wide$fasttt.lower_Black>
                                            ct.agency.wide$fasttt.upper_White,
                                          "1 White Drivers Disp. Neg.",
                                          "2 Approx. Equal"))
ct.agency.wide$dif.sr = ct.agency.wide$search_rate_white-
  ct.agency.wide$search_rate_black
ct.agency.wide$dif.sr.cat = ifelse(ct.agency.wide$dif.sr==0,
                                   "2 Apprrox. Equal",
                                   ifelse(ct.agency.wide$dif.sr>0,
                                          "1 White Drivers Disp. Neg.",
                                          "3 Black Drivers Disp. Neg."))
ct.agency.wide$dif.hr = ct.agency.wide$hit_rate_white-
  ct.agency.wide$hit_rate_black
ct.agency.wide$dif.hr.cat = ifelse(ct.agency.wide$dif.hr==0,
                                   "2 Approx. Equal",
                                   ifelse(ct.agency.wide$dif.hr<0,
                                          "1 White Drivers Disp. Neg.",
                                          "3 Black Drivers Disp. Neg."))
colnames(ct.agency.wide)[2:23] = paste(colnames(ct.agency.wide)[2:23],"noout",sep="_")

save(ct.agency.wide,file="Data/CT_WideAgencyData_50NoOut.RData")

ct.threshold.plot2 = ggplot(ct.agency.wide,aes(fasttt.mean_White_noout,fasttt.mean_Black_noout)) + 
  geom_point(size=1.5,shape=1) + theme_bw(base_size=12) +
  xlab(NULL) + ylab(NULL) +
  geom_abline(intercept = 0, slope = 1,color = "black",size = 0.5) +
  geom_abline(intercept = 0, slope = 2,color = "darkgray",linetype="dashed",size = 0.5) +
  geom_abline(intercept = 0, slope = 0.5,color = "darkgray",linetype="dashed",size = 0.5) +
  xlim(0,5)+ylim(0,5)+
  theme(panel.grid.major = element_line(colour = NA),
        panel.grid.minor = element_line(colour=NA),
        legend.position = "none") +
  ggtitle("Inferred Thresholds")
ct.threshold.plot2

#
# Combo
#

rm(list=ls())

# The combined measures
load("Data/CT_WideAgencyData_50Full.RData")
ct.combo = ct.agency.wide
load("Data/CT_WideAgencyData_50NoOut.RData")
ct.combo = cbind(ct.combo,ct.agency.wide[,2:23])

dim(ct.combo)
length(unique(apply(as.matrix(as.character(ct.combo$agency)),1,function(x){strsplit(x,"::",fixed=T)[[1]][1]})))
sum(ct.combo$stops_Black_noout)+sum(ct.combo$stops_White_noout)
sum(ct.combo$searches_Black_noout)+sum(ct.combo$searches_White_noout)

# Identify high disparity officers. 
load("Data/CT_WideOfficerData.RData")
ct.officer.wide$ratio = abs(ct.officer.wide$fasttt.mean_White/ct.officer.wide$fasttt.mean_Black)
summary(ct.officer.wide$ratio)
ct.officer.extreme.bw = boxplot(ct.officer.wide$ratio)
ct.officer.extreme.id = ct.officer.wide$officer[which(ct.officer.wide$ratio%in%
                                                        unique(ct.officer.extreme.bw$out))]

ct.agencies = unique(apply(as.matrix(as.character(ct.officer.wide$officer)),
                           1,
                           function(x){strsplit(x,"::",fixed=T)[[1]][1]}))
ct.agencies.extreme = unique(apply(as.matrix(as.character(ct.officer.wide$officer[which(ct.officer.wide$ratio%in%unique(ct.officer.extreme.bw$out))])),
                                   1,
                                   function(x){strsplit(x,"::",fixed=T)[[1]][1]}))

ct.combo$agency_any = ifelse(tolower(ct.combo$agency) %in% tolower(ct.agencies),1,0)
ct.combo$agency_out = ifelse(tolower(ct.combo$agency) %in% tolower(ct.agencies.extreme),1,0)

table(ct.combo$agency_any,ct.combo$agency_out)

ggplot(ct.combo %>% filter(agency_any==1),
       aes(fasttt.mean_White_full,fasttt.mean_Black_full,
           color=factor(agency_out),shape=factor(agency_out))) + 
  geom_point(size=1.5) + theme_bw(base_size=12) +
  xlab(NULL) + ylab(NULL) +
  geom_abline(intercept = 0, slope = 1,color = "black",size = 0.5) +
  geom_abline(intercept = 0, slope = 2,color = "darkgray",linetype="dashed",size = 0.5) +
  geom_abline(intercept = 0, slope = 0.5,color = "darkgray",linetype="dashed",size = 0.5) +
  xlim(0,5)+ylim(0,5)+
  theme(panel.grid.major = element_line(colour = NA),
        panel.grid.minor = element_line(colour=NA),
        legend.position = "none") +
  ggtitle("Inferred Thresholds")

ct.combo$fasttt.mean_Black_dif = ct.combo$fasttt.mean_Black_full - 
  ct.combo$fasttt.mean_Black_noout
ct.combo$fasttt.mean_White_dif = ct.combo$fasttt.mean_White_full - 
  ct.combo$fasttt.mean_White_noout
ct.combo$fasttt.black_dif_signif = ifelse(ct.combo$agency_out==1,
                                          ifelse(ct.combo$fasttt.mean_Black_dif<0,
                                                 ifelse(ct.combo$fasttt.mean_Black_full>
                                                          ct.combo$fasttt.lower_Black_noout & 
                                                          ct.combo$fasttt.mean_Black_noout<
                                                          ct.combo$fasttt.upper_White_full,1,0),
                                                 ifelse(ct.combo$fasttt.mean_Black_full<
                                                          ct.combo$fasttt.upper_Black_noout & 
                                                          ct.combo$fasttt.mean_Black_noout>
                                                          ct.combo$fasttt.lower_White_full,1,0)),
                                          NA)

ct.combo$thresh.black.signif = ifelse(abs(ct.combo$fasttt.mean_Black_dif) < 
                                        abs(ct.combo$fasttt.mean_Black_full-ct.combo$fasttt.upper_Black_full)&
                                        abs(ct.combo$fasttt.mean_Black_dif) < 
                                        abs(ct.combo$fasttt.mean_Black_noout-ct.combo$fasttt.upper_Black_noout),
                                      1,0)

ct.combo.sm = ct.combo %>% filter(agency_out==1)

mean(ct.combo.sm$fasttt.mean_White_full/ct.combo.sm$fasttt.mean_Black_full)
mean(ct.combo.sm$fasttt.mean_White_noout/ct.combo.sm$fasttt.mean_Black_noout)
t.test(ct.combo.sm$fasttt.mean_White_full/ct.combo.sm$fasttt.mean_Black_full,
       ct.combo.sm$fasttt.mean_White_noout/ct.combo.sm$fasttt.mean_Black_noout,paired = T)
sd(ct.combo.sm$fasttt.mean_White_full/ct.combo.sm$fasttt.mean_Black_full)

mean(ct.combo.sm$fasttt.mean_Black_full)
mean(ct.combo.sm$fasttt.mean_Black_noout)
t.test(ct.combo.sm$fasttt.mean_Black_full,
       ct.combo.sm$fasttt.mean_Black_noout,paired = T)
sd(ct.combo.sm$fasttt.mean_Black_full)

mean(ct.combo.sm$fasttt.mean_White_full)
mean(ct.combo.sm$fasttt.mean_White_noout)
t.test(ct.combo.sm$fasttt.mean_White_full,
       ct.combo.sm$fasttt.mean_White_noout,paired = T)
sd(ct.combo.sm$fasttt.mean_White_full)

t.test(ct.combo.sm$fasttt.mean_White_full-ct.combo.sm$fasttt.mean_White_noout,
       ct.combo.sm$fasttt.mean_Black_full-ct.combo.sm$fasttt.mean_Black_noout,
       paired = T)

prop.table(table(ct.combo$dif.sr.cat_noout))
prop.table(table(ct.combo$dif.tt.cat_noout))

sd(ct.combo.sm$fasttt.mean_White_full)
sd(ct.combo.sm$fasttt.mean_Black_full)
sd(ct.combo.sm$fasttt.mean_White_noout/ct.combo.sm$fasttt.mean_Black_noout)

ct_noout_sr = ggplot(ct.combo,aes(search_rate_white_noout,
                                  search_rate_black_noout)) + 
  geom_point(size=1.5) + theme_bw(base_size=12) +
  xlab(NULL) + ylab(NULL) +
  geom_abline(intercept = 0, slope = 1,
              color = "darkgray",size = 0.5) +
  geom_abline(intercept = 0, slope = 2,
              color = "darkgray",linetype="dashed",size = 0.5) +
  geom_abline(intercept = 0, slope = 0.5,
              color = "darkgray",linetype="dashed",size = 0.5) +
  xlim(0,0.25)+ylim(0,0.25)+
  theme(panel.grid.major = element_line(colour = NA),
        panel.grid.minor = element_line(colour=NA),
        legend.position = "none")
ggsave(plot=ct_noout_sr,filename = "Figures/Figure2a.eps")
png("Figures/CT_SR_NoOut.png",547,441)
ct_noout_sr
dev.off()

ct_noout_tt = ggplot(ct.combo,aes(fasttt.mean_White_noout,
                                  fasttt.mean_Black_noout)) + 
  geom_point(size=1.5) + theme_bw(base_size=12) +
  xlab(NULL) + ylab(NULL) +
  geom_abline(intercept = 0, slope = 1,
              color = "darkgray",size = 0.5) +
  geom_abline(intercept = 0, slope = 2,
              color = "darkgray",linetype="dashed",size = 0.5) +
  geom_abline(intercept = 0, slope = 0.5,
              color = "darkgray",linetype="dashed",size = 0.5) +
  xlim(0,4)+ylim(0,4)+
  theme(panel.grid.major = element_line(colour = NA),
        panel.grid.minor = element_line(colour=NA),
        legend.position = "none")
ggsave(plot=ct_noout_tt,filename = "Figures/Figure2b.eps")
png("Figures/CT_TT_NoOut.png",547,441)
ct_noout_tt
dev.off()

